* GECM w/ multiple predictors DGP

// Basics
clear all
set matsize 11000
 
parallel setclusters 8

capture program drop ecm
program define ecm, rclass

// Monte Carlo Basics
version 12.1
syntax [, obs(integer 1) rho_e(real 0.5)]

drop _all

*args obs phi_u beta_x

set obs `obs'

// TS Basics
egen time = fill(0 1 2)
tsset time

// Data Generating Process
gen x1 = rnormal() if time == 0
replace x1 = l.x1 + rnormal() if time > 0

gen x2 = rnormal() if time == 0
replace x2 = l.x2 + rnormal() if time > 0

scalar beta1 = 2
scalar gamma1 = 1.5
scalar beta2 = 1
scalar gamma2 = 1.25

gen y = rnormal()
replace y = `rho_e'*l.y + beta1*d.x1 + gamma1*l.x1 + beta2*d.x2 + gamma2*l.x2 + rnormal() if time > 0

dfuller y
return scalar stat_y = r(p) < 0.05 

// Model estimation
local gecm_true_phi = `rho_e' - 1
// Model Estimation
* ADL (1,1)
reg y L.y x1 L.x1 x2 L.x2
return scalar phi_adl = _b[L1.y]
scalar phi_adl_test_t1 = (_b[L1.y] - `rho_e') / _se[L1.y]
return scalar phi_adl_t1 = 2*ttail(e(df_r),abs(phi_adl_test_t1)) < 0.05
scalar phi_adl_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_adl_power = 2*ttail(e(df_r),abs(phi_adl_test_power)) < 0.05
return scalar b1_adl = _b[x1]
scalar b1_adl_test_t1 = (_b[x1] - beta1) / _se[x1]
return scalar b1_adl_t1 = 2*ttail(e(df_r),abs(b1_adl_test_t1)) < 0.05
scalar b1_adl_test_power = (_b[x1] - 0) / _se[x1]
return scalar b1_adl_power = 2*ttail(e(df_r),abs(b1_adl_test_power)) < 0.05
return scalar g1_adl = _b[L1.x1]
scalar g1_adl_test_t1 = (_b[L1.x1] - (gamma1-beta1)) / _se[L1.x1]
return scalar g1_adl_t1 = 2*ttail(e(df_r),abs(g1_adl_test_t1)) < 0.05
scalar g1_adl_test_power = (_b[L1.x1] - 0) / _se[L1.x1]
return scalar g1_adl_power = 2*ttail(e(df_r),abs(g1_adl_test_power)) < 0.05
return scalar b2_adl = _b[x2]
scalar b2_adl_test_t1 = (_b[x2] - beta2) / _se[x2]
return scalar b2_adl_t1 = 2*ttail(e(df_r),abs(b2_adl_test_t1)) < 0.05
scalar b2_adl_test_power = (_b[x2] - 0) / _se[x2]
return scalar b2_adl_power = 2*ttail(e(df_r),abs(b2_adl_test_power)) < 0.05
return scalar g2_adl = _b[L1.x2]
scalar g2_adl_test_t1 = (_b[L1.x2] - (gamma2 - beta2)) / _se[L1.x2]
return scalar g2_adl_t1 = 2*ttail(e(df_r),abs(g2_adl_test_t1)) < 0.05
scalar g2_adl_test_power = (_b[L1.x2] - 0) / _se[L1.x2]
return scalar g2_adl_power = 2*ttail(e(df_r),abs(g2_adl_test_power)) < 0.05
return scalar lrm_adl_x1 = (_b[x1] + _b[L1.x1]) /(1-_b[L1.y])
nlcom ((_b[x1] + _b[L1.x1]) /(1-_b[L1.y])) * (1e+3)
matrix define var_lrm_adl_x1 = r(V)
scalar var_lrm_adl_x1_se  = sqrt(var_lrm_adl_x1[1,1]) / (1e+3)
scalar lrm_adl_x1_test_t1 = (((_b[x1] + _b[L1.x1]) /(1-_b[L1.y])) - ((beta1 + (gamma1 - beta1))/(1-`rho_e'))) /var_lrm_adl_x1_se
return scalar lrm_adl_x1_t1 = 2*ttail(e(df_r),abs(lrm_adl_x1_test_t1)) < 0.05
scalar lrm_adl_x1_test_power = (((_b[x1] + _b[L1.x1]) /(1-_b[L1.y])) - 0) /var_lrm_adl_x1_se
return scalar lrm_adl_x1_power = 2*ttail(e(df_r),abs(lrm_adl_x1_test_power)) < 0.05
return scalar lrm_adl_x2 = (_b[x2] + _b[L1.x2]) /(1-_b[L1.y])
nlcom ((_b[x2] + _b[L1.x2]) /(1-_b[L1.y])) * (1e+3)
matrix define var_lrm_adl_x2 = r(V)
scalar var_lrm_adl_x2_se  = sqrt(var_lrm_adl_x2[1,1]) / (1e+3)
scalar lrm_adl_x2_test_t1 = (((_b[x2] + _b[L1.x2]) /(1-_b[L1.y])) - ((- gamma2)/(`gecm_true_phi'))) /var_lrm_adl_x2_se
return scalar lrm_adl_x2_t1 = 2*ttail(e(df_r),abs(lrm_adl_x2_test_t1)) < 0.05
scalar lrm_adl_x2_test_power = (((_b[x2] + _b[L1.x2]) /(1-_b[L1.y])) - 0) /var_lrm_adl_x2_se
return scalar lrm_adl_x2_power = 2*ttail(e(df_r),abs(lrm_adl_x2_test_power)) < 0.05

* GECM
reg D.y L.y D.x1 L.x1 D.x2 L.x2
return scalar phi_gecm = _b[L1.y]
local gecm_true_phi = `rho_e' - 1
scalar phi_gecm_test_t1 = (_b[L1.y] - `gecm_true_phi') / _se[L1.y]
return scalar phi_gecm_t1 = 2*ttail(e(df_r),abs(phi_gecm_test_t1)) < 0.05
scalar phi_gecm_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_gecm_power = 2*ttail(e(df_r),abs(phi_gecm_test_power)) < 0.05
return scalar b1_gecm  = _b[D1.x1]
scalar b1_gecm_test_t1= (_b[D1.x1] - beta1) / _se[D1.x1]
return scalar b1_gecm_t1 = 2*ttail(e(df_r),abs(b1_gecm_test_t1)) < 0.05
scalar b1_gecm_test_power = (_b[D1.x1] - 0) / _se[D1.x1]
return scalar b1_gecm_power = 2*ttail(e(df_r),abs(b1_gecm_test_power)) < 0.05
return scalar g1_gecm  = _b[L1.x1]
scalar g1_gecm_test_t1 = (_b[L1.x1] - (gamma1)) / _se[L1.x1]
return scalar g1_gecm_t1 = 2*ttail(e(df_r),abs(g1_gecm_test_t1)) < 0.05
scalar g1_gecm_test_power = (_b[L1.x1] - 0) / _se[L1.x1]
return scalar g1_gecm_power = 2*ttail(e(df_r),abs(g1_gecm_test_power)) < 0.05
return scalar b2_gecm  = _b[D1.x2]
scalar b2_gecm_test_t1= (_b[D1.x2] - beta2) / _se[D1.x2]
return scalar b2_gecm_t1 = 2*ttail(e(df_r),abs(b2_gecm_test_t1)) < 0.05
scalar b2_gecm_test_power = (_b[D1.x2] - 0) / _se[D1.x2]
return scalar b2_gecm_power = 2*ttail(e(df_r),abs(b2_gecm_test_power)) < 0.05
return scalar g2_gecm  = _b[L1.x2]
scalar g2_gecm_test_t1 = (_b[L1.x2] - ( gamma2)) / _se[L1.x2]
return scalar g2_gecm_t1 = 2*ttail(e(df_r),abs(g2_gecm_test_t1)) < 0.05
scalar g2_gecm_test_power = (_b[L1.x2] - 0) / _se[L1.x2]
return scalar g2_gecm_power = 2*ttail(e(df_r),abs(g2_gecm_test_power)) < 0.05
return scalar lrm_gecm_x1 = -_b[L1.x1] / _b[L1.y]
nlcom (-_b[L1.x1] / _b[L1.y])* (1e+3)
matrix define var_lrm_gecm_x1 = r(V)
scalar var_lrm_gecm_x1_se  = sqrt(var_lrm_gecm_x1[1,1]) / (1e+3)
scalar lrm_gecm_x1_test_t1 = ((-_b[L1.x1] / _b[L1.y]) - ((- gamma1)/(`gecm_true_phi'))) /var_lrm_gecm_x1_se
return scalar lrm_gecm_x1_t1 = 2*ttail(e(df_r),abs(lrm_gecm_x1_test_t1)) < 0.05
scalar lrm_gecm_x1_test_power = ((-_b[L1.x1] / _b[L1.y]) - 0) /var_lrm_gecm_x1_se
return scalar lrm_gecm_x1_power = 2*ttail(e(df_r),abs(lrm_gecm_x1_test_power)) < 0.05
return scalar lrm_gecm_x2 = -_b[L1.x2] / _b[L1.y]
nlcom (-_b[L1.x2] / _b[L1.y])* (1e+3)
matrix define var_lrm_gecm_x2 = r(V)
scalar var_lrm_gecm_x2_se  = sqrt(var_lrm_gecm_x2[1,1]) / (1e+3)
scalar lrm_gecm_x2_test_t1 = ((-_b[L1.x2] / _b[L1.y]) - ((- gamma2)/(`gecm_true_phi'))) /var_lrm_gecm_x2_se
return scalar lrm_gecm_x2_t1 = 2*ttail(e(df_r),abs(lrm_gecm_x2_test_t1)) < 0.05
scalar lrm_gecm_x2_test_power = ((-_b[L1.x2] / _b[L1.y]) - 0) /var_lrm_gecm_x2_se
return scalar lrm_gecm_x2_power = 2*ttail(e(df_r),abs(lrm_gecm_x2_test_power)) < 0.05


* ADL (2,2)
reg y L.y L2.y x1 L.x1 L2.x1 x2 L.x2 L2.x2
return scalar phi1_adl22 = _b[L1.y]
scalar phi1_adl22_test_t1 = (_b[L1.y] - `rho_e') / _se[L1.y]
return scalar phi1_adl22_t1 = 2*ttail(e(df_r),abs(phi1_adl22_test_t1)) < 0.05
scalar phi1_adl22_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi1_adl22_power = 2*ttail(e(df_r),abs(phi1_adl22_test_power)) < 0.05
return scalar phi2_adl22 = _b[L2.y]
scalar phi2_adl22_test = (_b[L2.y] - 0) / _se[L2.y]
return scalar phi2_adl22_t1 = 2*ttail(e(df_r),abs(phi2_adl22_test)) < 0.05
return scalar b1_adl22 = _b[x1]
scalar b1_adl22_test_t1 = (_b[x1] - beta1) / _se[x1]
return scalar b1_adl22_t1 = 2*ttail(e(df_r),abs(b1_adl22_test_t1)) < 0.05
scalar b1_adl22_test_power = (_b[x1] - 0) / _se[x1]
return scalar b1_adl22_power = 2*ttail(e(df_r),abs(b1_adl22_test_power)) < 0.05
return scalar g1_adl22 = _b[L1.x1]
scalar g1_adl22_test_t1 = (_b[L1.x1] - (gamma1 -beta1)) / _se[L1.x1]
return scalar g1_adl22_t1 = 2*ttail(e(df_r),abs(g1_adl22_test_t1)) < 0.05
scalar g1_adl22_test_power = (_b[L1.x1] - 0) / _se[L1.x1]
return scalar g1_adl22_power = 2*ttail(e(df_r),abs(g1_adl22_test_power)) < 0.05
return scalar b2_adl22 = _b[x2]
scalar b2_adl22_test_t1 = (_b[x2] - beta2) / _se[x2]
return scalar b2_adl22_t1 = 2*ttail(e(df_r),abs(b2_adl22_test_t1)) < 0.05
scalar b2_adl22_test_power = (_b[x2] - 0) / _se[x2]
return scalar b2_adl22_power = 2*ttail(e(df_r),abs(b2_adl22_test_power)) < 0.05
return scalar g2_adl22 = _b[L1.x2]
scalar g2_adl22_test_t1 = (_b[L1.x2] - (gamma2-beta2)) / _se[L1.x2]
return scalar g2_adl22_t1 = 2*ttail(e(df_r),abs(g2_adl22_test_t1)) < 0.05
scalar g2_adl22_test_power = (_b[L1.x2] - 0) / _se[L1.x2]
return scalar g2_adl22_power = 2*ttail(e(df_r),abs(g2_adl22_test_power)) < 0.05
return scalar lrm_adl22_x1 = (_b[x1] + _b[L1.x1] + _b[L2.x1]) /(1-_b[L1.y] - _b[L2.y])
nlcom ((_b[x1] + _b[L1.x1] + _b[L2.x1]) /(1-_b[L1.y]- _b[L2.y])) * (1e+3)
matrix define var_lrm_adl22_x1 = r(V)
scalar var_lrm_adl22_x1_se  = sqrt(var_lrm_adl22_x1[1,1]) / (1e+3)
scalar lrm_adl22_x1_test_t1 = (((_b[x1] + _b[L1.x1] + _b[L2.x1]) /(1-_b[L1.y] - _b[L2.y])) - ((beta1 + (gamma1 - beta1))/(1-`rho_e'))) /var_lrm_adl22_x1_se
return scalar lrm_adl22_x1_t1 = 2*ttail(e(df_r),abs(lrm_adl22_x1_test_t1)) < 0.05
scalar lrm_adl22_x1_test_power = (((_b[x1] + _b[L1.x1]+ _b[L2.x1]) /(1-_b[L1.y]-_b[L2.y])) - 0) /var_lrm_adl22_x1_se
return scalar lrm_adl22_x1_power = 2*ttail(e(df_r),abs(lrm_adl22_x1_test_power)) < 0.05
return scalar lrm_adl22_x2 = (_b[x2] + _b[L1.x2] + _b[L2.x2]) /(1-_b[L1.y]-_b[L2.y])
nlcom ((_b[x2] + _b[L1.x2] + _b[L2.x2]) /(1-_b[L1.y]-_b[L2.y])) * (1e+3)
matrix define var_lrm_adl22_x2 = r(V)
scalar var_lrm_adl22_x2_se  = sqrt(var_lrm_adl22_x2[1,1]) / (1e+3)
scalar lrm_adl22_x2_test_t1 = (((_b[x2] + _b[L1.x2] + _b[L2.x2]) /(1-_b[L1.y]-_b[L2.y])) - (((beta2 + (gamma2 - beta2))/(1-`rho_e')))) /var_lrm_adl22_x2_se
return scalar lrm_adl22_x2_t1 = 2*ttail(e(df_r),abs(lrm_adl22_x2_test_t1)) < 0.05
scalar lrm_adl22_x2_test_power = (((_b[x2] + _b[L1.x2]+ _b[L2.x2]) /(1-_b[L1.y]-_b[L2.y])) - 0) /var_lrm_adl22_x2_se
return scalar lrm_adl22_x2_power = 2*ttail(e(df_r),abs(lrm_adl22_x2_test_power)) < 0.05
return scalar b3_adl22 = _b[L2.x1]
scalar b3_adl22_test = (_b[L2.x1] - 0) / _se[L2.x1]
return scalar b3_adl22_t1 = 2*ttail(e(df_r),abs(b3_adl22_test)) < 0.05
return scalar b4_adl22 = _b[L2.x2]
scalar b4_adl22_test = (_b[L2.x2] - 0) / _se[L2.x2]
return scalar b4_adl22_t1 = 2*ttail(e(df_r),abs(b4_adl22_test)) < 0.05

* ADL (3,3)
reg y L.y L2.y L3.y x1 L.x1 L2.x1 L3.x1 x2 L.x2 L2.x2 L3.x2
return scalar phi1_adl33 = _b[L1.y]
scalar phi1_adl33_test_t1 = (_b[L1.y] - `rho_e') / _se[L1.y]
return scalar phi1_adl33_t1 = 2*ttail(e(df_r),abs(phi1_adl33_test_t1)) < 0.05
scalar phi1_adl33_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi1_adl33_power = 2*ttail(e(df_r),abs(phi1_adl33_test_power)) < 0.05
return scalar phi2_adl33 = _b[L2.y]
scalar phi2_adl33_test = (_b[L2.y] - 0) / _se[L2.y]
return scalar phi2_adl33_t1 = 2*ttail(e(df_r),abs(phi2_adl33_test)) < 0.05
return scalar phi3_adl33 = _b[L3.y]
scalar phi3_adl33_test = (_b[L3.y] - 0) / _se[L3.y]
return scalar phi3_adl33_t1 = 2*ttail(e(df_r),abs(phi3_adl33_test)) < 0.05
return scalar b1_adl33 = _b[x1]
scalar b1_adl33_test_t1 = (_b[x1] - beta1) / _se[x1]
return scalar b1_adl33_t1 = 2*ttail(e(df_r),abs(b1_adl33_test_t1)) < 0.05
scalar b1_adl33_test_power = (_b[x1] - 0) / _se[x1]
return scalar b1_adl33_power = 2*ttail(e(df_r),abs(b1_adl33_test_power)) < 0.05
return scalar g1_adl33 = _b[L1.x1]
scalar g1_adl33_test_t1 = (_b[L1.x1] - (gamma1-beta1)) / _se[L1.x1]
return scalar g1_adl33_t1 = 2*ttail(e(df_r),abs(g1_adl33_test_t1)) < 0.05
scalar g1_adl33_test_power = (_b[L1.x1] - 0) / _se[L1.x1]
return scalar g1_adl33_power = 2*ttail(e(df_r),abs(g1_adl33_test_power)) < 0.05
return scalar b2_adl33 = _b[x2]
scalar b2_adl33_test_t1 = (_b[x2] - beta2) / _se[x2]
return scalar b2_adl33_t1 = 2*ttail(e(df_r),abs(b2_adl33_test_t1)) < 0.05
scalar b2_adl33_test_power = (_b[x2] - 0) / _se[x2]
return scalar b2_adl33_power = 2*ttail(e(df_r),abs(b2_adl33_test_power)) < 0.05
return scalar g2_adl33= _b[L1.x2]
scalar g2_adl33_test_t1 = (_b[L1.x2] - (gamma2-beta2)) / _se[L1.x2]
return scalar g2_adl33_t1 = 2*ttail(e(df_r),abs(g2_adl33_test_t1)) < 0.05
scalar g2_adl33_test_power = (_b[L1.x2] - 0) / _se[L1.x2]
return scalar g2_adl33_power = 2*ttail(e(df_r),abs(g2_adl33_test_power)) < 0.05
return scalar lrm_adl33_x1 = (_b[x1] + _b[L1.x1] + _b[L2.x1] + _b[L3.x1]) /(1-_b[L1.y] - _b[L2.y] - _b[L3.y])
nlcom ((_b[x1] + _b[L1.x1] + _b[L2.x1] + _b[L3.x1]) /(1-_b[L1.y] - _b[L2.y] - _b[L3.y])) * (1e+3)
matrix define var_lrm_adl33_x1 = r(V)
scalar var_lrm_adl33_x1_se  = sqrt(var_lrm_adl33_x1[1,1]) / (1e+3)
scalar lrm_adl33_x1_test_t1 = (((_b[x1] + _b[L1.x1] + _b[L2.x1] + _b[L3.x1]) /(1-_b[L1.y] - _b[L2.y]- _b[L3.y])) - ((beta1 + (gamma1-beta1))/(1-`rho_e'))) /var_lrm_adl22_x1_se
return scalar lrm_adl33_x1_t1 = 2*ttail(e(df_r),abs(lrm_adl33_x1_test_t1)) < 0.05
scalar lrm_adl33_x1_test_power = (((_b[x1] + _b[L1.x1]+ _b[L2.x1] + _b[L3.x1]) /(1-_b[L1.y]-_b[L2.y]- _b[L3.y])) - 0) /var_lrm_adl33_x1_se
return scalar lrm_adl33_x1_power = 2*ttail(e(df_r),abs(lrm_adl33_x1_test_power)) < 0.05
return scalar lrm_adl33_x2 = (_b[x2] + _b[L1.x2] + _b[L2.x2] + _b[L3.x2]) /(1-_b[L1.y]-_b[L2.y]- _b[L3.y])
nlcom ((_b[x2] + _b[L1.x2] + _b[L2.x2] + _b[L3.x2]) /(1-_b[L1.y]-_b[L2.y]- _b[L3.y])) * (1e+3)
matrix define var_lrm_adl33_x2 = r(V)
scalar var_lrm_adl33_x2_se  = sqrt(var_lrm_adl33_x2[1,1]) / (1e+3)
scalar lrm_adl33_x2_test_t1 = (((_b[x2] + _b[L1.x2]+ _b[L2.x2] + _b[L3.x2]) /(1-_b[L1.y]-_b[L2.y]- _b[L3.y])) - ((beta2 + (gamma2-beta2))/(1-`rho_e'))) /var_lrm_adl33_x2_se
return scalar lrm_adl33_x2_t1 = 2*ttail(e(df_r),abs(lrm_adl33_x2_test_t1)) < 0.05
scalar lrm_adl33_x2_test_power = (((_b[x2] + _b[L1.x2]+ _b[L2.x2] + _b[L3.x2]) /(1-_b[L1.y]-_b[L2.y]- _b[L3.y])) - 0) /var_lrm_adl33_x2_se
return scalar lrm_adl33_x2_power = 2*ttail(e(df_r),abs(lrm_adl33_x2_test_power)) < 0.05
return scalar b3_adl33 = _b[L2.x1]
scalar b3_adl33_test = (_b[L2.x1] - 0) / _se[L2.x1]
return scalar b3_adl33_t1 = 2*ttail(e(df_r),abs(b3_adl33_test)) < 0.05
return scalar b4_adl33 = _b[L2.x2]
scalar b4_adl33_test = (_b[L2.x2] - 0) / _se[L2.x2]
return scalar b4_adl33_t1 = 2*ttail(e(df_r),abs(b4_adl33_test)) < 0.05
return scalar b5_adl33 = _b[L3.x1]
scalar b5_adl33_test = (_b[L3.x1] - 0) / _se[L3.x1]
return scalar b5_adl33_t1 = 2*ttail(e(df_r),abs(b5_adl33_test)) < 0.05
return scalar b6_adl33 = _b[L3.x2]
scalar b6_adl33_test = (_b[L3.x2] - 0) / _se[L3.x2]
return scalar b6_adl33_t1 = 2*ttail(e(df_r),abs(b6_adl33_test)) < 0.05

end
